clear;
load('options_BE.mat')
load('Data_FR2_r.mat', 'Data_FR2_r')
nsect = length(Data_FR2_r);
%nsect=1
kjInd_FR2 = zeros(nsect,1);

for s = 1:nsect
    sigma=Data_FR2_r(s,3);    
    %sigma=1.0
    k_cool = @(k_j)sd_loge_int(k_j, sigma);
    x = fsolve(k_cool,5, options_BE)
    kjInd_FR2(s,1)=x;
    save('kjInd_FR2.mat', 'kjInd_FR2');
end



save('kjInd_FR2.mat', 'kjInd_FR2');


kappa1 = zeros(nsect,1);
kappa2 = zeros(nsect,1);
theta = zeros(nsect,1);

for s = 1:nsect
fun3 = @(z, k) (1-(z.^2)).*(z.*(exp(z))).^k .*exp(z);
Integral3 = integral(@(z)fun3(z, kjInd_FR2(s,1)),0,1);

kappa1(s,1)=kjInd_FR2(s,1).*exp(-kjInd_FR2(s,1)-1).*Integral3;
kappa2(s,1)=kappa1(s,1)./kjInd_FR2(s,1);
theta(s,1)=(kappa2(s,1).*(kjInd_FR2(s,1)+1))./((1./(kjInd_FR2(s,1)+1))-(kappa1(s,1)+kappa2(s,1))  );
end
thetajInd_FR2=theta;
save('thetajInd_FR2.mat', 'thetajInd_FR2');


init = ones(nsect,1);
CC=Data_FR2_r(:,4);
AA=theta;
ns=s;
beta_cool= @(B)eqsystem(AA, B, CC, ns);
    x = fsolve(beta_cool, init, options_BE)
    betajInd_FR2=x;
    save('betajInd_FR2.mat', 'betajInd_FR2');
    
    
intra_dist_FR2 = zeros(nsect,1);
intra_dist_FR2=100.*( ((kappa2.*((kjInd_FR2+1).^2)).^(-1./(kjInd_FR2+1))) -1  );
trick=betajInd_FR2'*thetajInd_FR2;
inter_dist_FR2 = zeros(nsect,1);
inter_dist_FR2=100.*((thetajInd_FR2./(trick))-1);

        save('inter_dist_FR2.mat', 'inter_dist_FR2');
        save('intra_dist_FR2.mat', 'intra_dist_FR2');